Model reduction =============== A heterogeneous-agent model written with ``@heterogeneity_axis`` expands into one copy of every individual variable per grid cell, plus the distribution and the interpolation weights: a household block of a few hundred cells already gives a model of thousands of equations, and every solve and every evaluation of the likelihood pays for all of them. The option ``solve_reduction`` solves and filters a much smaller system instead, built once from the model. .. code-block:: matlab R = rise.engine.dsge_tools.reduction.build(m); % once, at reference parameters m = set(m, 'solve_reduction', R); % then as usual [~, LL] = filter(m, 'data', db); m = estimate(m, ...); The option is off by default (``[]``). Setting it back to ``[]`` restores the full solve. What the reduction does ----------------------- The variables are split in two. The **kept** variables -- the aggregates, the observables, and anything that drives a time-varying transition probability -- stay exactly as they are. The **reduced** variables -- by default everything the heterogeneity expansion generated: policies, multipliers, the distribution, the interpolation weights -- are written as their steady state plus a basis times a few coordinates, .. math:: x_h = \bar x_h + V z . :math:`V` is block diagonal: one orthonormal block per family of variables (consumption, savings, the distribution, ...) and per RISE type (static, predetermined, forward-looking), so each coordinate has the type of its block and the reduced system keeps the timing structure of the full one. With :math:`x = E x_r`, :math:`E = \mathrm{blkdiag}(I, V)`, the equations are projected with .. math:: W' = E' A_0^{-1}, where :math:`A_0` is the matrix that maps equation residuals into the impact responses of the variables at the reference solution. The projected equations are the full equations expressed in the units of the variables and restricted to the span of the basis. If the solution of the full model lies in that span, every one of its paths solves the reduced system: the reduction is then **exact, whatever the parameters of the equations**. Its accuracy is a property of the basis only. The reduced problem is solved by the usual solvers, filtered in its own (small) state space, and its solution is re-inflated to the full variables, so ``irf``, ``simulate`` and the smoothed variables report every variable of the model. The basis --------- ``build`` solves the model at its current parameters (first order) and computes, group by group, a proper orthogonal decomposition of * the covariance of the variables under each shock alone, in each regime, normalized to unit trace (so the scale of the shocks does not matter); * the responses to an anticipated forcing :math:`u_t=\rho u_{t-1}+e_t` of every equation that involves no reduced variable, for persistences up to 0.999. The model's own shocks only excite the directions the aggregate dynamics produce at the reference parameters; the forcing spans the responses of the household block to the price paths other parameters produce. Without it the basis fails as soon as persistences or policy coefficients move; * at ``solve_order`` 2 and above, the responses to the second-order forcing terms (see below). The dominant directions of each source are selected separately (option ``tol`` for the shocks, ``forcing_tol`` for the forcing) and merged. Tighter is not always better: directions far down the spectrum are numerically noisy. The basis is computed **once**. It is then used for any other parameter vector, which is what makes the reduction useful in estimation: per evaluation only the projection and the small solve are done. It must be built at the steady state it will be used at: the reduced coordinates are deviations from it, and a solve at a steady state where the reduced variables have moved is refused. Regime switching ---------------- One basis serves all regimes -- it spans the covariances of every regime -- and the same projection is applied to every pair of regimes. The transition probabilities that the parser multiplies into each block pass through it unchanged, so constant and time-varying transition probabilities go through exactly the same path. A variable that enters a time-varying transition probability is kept; the filters evaluate the probabilities on the full vector implied by the reduced state. Mean-square stability is checked on the reduced solution: the full solution lives in the span of the basis, so the second-moment operator of the full states has the same nonzero eigenvalues as that of the reduced states, which is a fraction of its size. Occasionally-binding constraints -------------------------------- Constraints are enforced by the filters exactly as for the full model, on the individual variables themselves: the constraint-enforcing forecaster of the full model is evaluated on the vector rebuilt from the reduced state (steady state plus basis times coordinates), and the enforcing shocks go through the projection like any other shock. The enforcers need a decision. There is one per constraint -- one per cell for a borrowing limit -- and each moves its own cell: putting all their responses in the basis leaves the household block unreduced. The option ``enforcers`` of ``build`` chooses which ones the basis represents: ``'auto'`` (default) those whose constraint switches status in ``enforcer_draws`` (5000) draws from the reference ergodic distribution of the variables, the shocks scaled by ``enforcer_scale`` (2): the constraints within about eight standard deviations of switching; ``'all'`` every one; ``'none'`` none. Whatever the choice, every constraint is enforced; the option only decides which enforcing responses the reduced system represents exactly. A constraint whose enforcing response the basis leaves out is enforced through the part of the response that the basis spans, and the filters, which update through the constrained step (see :doc:`/WorkingWithAModel/Filtering`), need that part to move the constraint. Under ``'A0'`` it can barely move it: on the model of the reduction tests (an Aiyagari economy with a borrowing limit and volatility regimes), a DCT and copula basis with ``'none'`` amplifies the slopes of the step up to :math:`2\times 10^{19}` times where the limit switches (at most 1.09 in the full model), the filters leave such steps to the Kalman update with a warning (``rise:filter:illConditionedConstrainedUpdate``) and mark them in their output, and the log-likelihood is -54.9 against -17.8 for the full model; with the default selection it is the full model's to :math:`10^{-4}`. ``'none'`` is safe only when no constraint binds in the sample; otherwise keep the default under ``'A0'``. The automatic selection is checked. If the reduced model built with it has no unique stable solution at the reference (its root count, ``R.roots``), the responses it left out are directions the solution needs: ``build`` rebuilds the reduction with every enforcer in the basis, says so in a warning (``rise:reduction:enforcersAll``), and records it in ``R.enforcers.mode``. On the one-asset New Keynesian HANK with the zero lower bound of the HANK course (12 x 3 grid, 37 constraints), the automatic selection keeps 6 enforcers and gives a reduced model with too few stable roots; with all 37 the reduced log-likelihood equals the full one to :math:`4\times 10^{-8}`. More generally, any reduction whose reduced model is not determinate at the reference is reported when it is built (``rise:reduction:notDeterminate``), before a solve or a filter fails on it. Higher orders ------------- The substitution :math:`x = \bar x + E x_r` and the projection are linear, so the derivatives of order :math:`k` of the reduced model are those of the full model contracted with :math:`W'` on the equations and with the basis on each of their :math:`k` variable indices. The reduced problem is solved by the ``hops`` engine (orders 1 to 5) like any model, and every order of the solution is re-inflated. The basis is built for the order the model is set to solve: set ``solve_order`` before calling ``build``. At order 2 and above, the quadratic terms of the solution are the responses of the first-order operator to the forcing :math:`f_{vv}(a\otimes b)` (pairs of first-order responses, and the precautionary term); their responses, for the stochastic shocks and the dominant state directions (``second_order_directions``), computed from the full model's second derivatives at the reference steady state, are added to the decomposition. This is the enrichment of Reiter (2023), *State Reduction and Second-order Perturbations of Heterogeneous Agent Models*. The full model is never solved at a higher order, but its second derivatives are needed: build it with ``'max_deriv_order', 2`` (or more). A reduction built for a lower order than the one solved is refused. Higher-order terms grow with a power of the size of the state: a reduction that is small enough at first order can still be too large at second or third order. The reduced derivatives of order :math:`k` are contracted from the nonzero entries of the full ones, by blocks of equations, with :math:`W'` applied block by block, so the only dense array is the reduced one: the reduced equations times the reduced differentiation variables (leads, current values, lags, shocks) to the power :math:`k`. The full model's size does not enter. Beyond the option ``max_dense`` of ``build`` (:math:`3\times 10^8` elements by default) the contraction refuses rather than exhaust memory: use a smaller basis (a larger ``tol``) or, if the machine has the memory, raise ``max_dense``. Supplied bases and projections ------------------------------ The decomposition above needs the solution of the full model at the reference parameters, and it adapts the basis to that solution. Two options of ``build`` give the other route, where the basis is fixed in advance and the full model is never solved: ``basis`` one field per family, holding its basis: a matrix with one row per member of the family (in the order of the grid cells, axis 1 fastest) and one column per coordinate. It replaces the decomposition for that family. ``'exact'`` (or ``[]``) keeps the family whole. ``projection`` how the equations are reduced. ``'A0'`` (the default) is the projection described above. ``'galerkin'`` pairs every variable with an equation and projects the equations of a family on the family's own basis. ``'equation_basis'`` does the same with a basis supplied for the equations (option ``equation_basis``). With ``'galerkin'`` and ``'equation_basis'`` only the steady state and the derivatives of the full model are computed. The reduction then applies to a full model that is too large to solve, or whose linearization has no unique stable solution. Two helpers build parameter-free bases. ``dct_basis(sizes, k)`` The ``k`` tensor DCT-II vectors of lowest frequency on a grid. With ``'energy', tol, 'reference', X``, it keeps the coefficients that carry a share ``1-tol`` of the energy of the columns of ``X``: for instance the steady-state policy or value function, as in Bayer and Luetticke (2020). ``marginal_copula_basis(pmf)`` The linearized map from perturbations of the marginal distributions, and optionally a few perturbations of the copula, to the joint distribution. The copula is held at its steady-state value. Every column sums to zero, so the mass never moves; a marginal perturbation leaves the other marginals unchanged. The basis is not orthonormal: the reduction maps a distribution back to its coordinates with the pseudo-inverse. A Bayer-Luetticke-type reduction of a one-asset model: .. code-block:: matlab het = m.model_data.heterogeneity; % grid, families, names D = cellfun(@(n) ss.(n), het.dist_name); % steady-state distribution B = struct(); B.c = rise.engine.dsge_tools.reduction.dct_basis(het.sizes, 8); B.D = rise.engine.dsge_tools.reduction.marginal_copula_basis( ... reshape(D, het.sizes), 'copula', 5); B.ap = 'exact'; B.mu = 'exact'; B.theta_a = 'exact'; R = rise.engine.dsge_tools.reduction.build(m, 'basis', B, ... 'projection', 'galerkin'); R.roots % root count of the reduced model mr = set(m, 'solve_reduction', R); ``R.roots`` reports, for each regime, the number of stable roots of the reduced first-order system against the number of predetermined variables, with a verdict: ``determinate``, ``indeterminate`` or ``no stable solution``. **Pairing equations with variables.** A Galerkin projection needs to know which equations belong to a family. By default every variable is paired with an equation by a matching that maximizes the product of the current-period derivatives. A cell's law of motion is then paired with its own mass, and the mass equation with one of the cells. Option ``equation_families`` imposes the pairs of a family, for instance the Euler equations for a value-function family: .. code-block:: matlab eqs = get(m, 'equations'); R = rise.engine.dsge_tools.reduction.build(m, 'basis', B, ... 'projection', 'galerkin', 'equation_families', struct('vm', euler_rows)); When a basis has fewer columns than rows, the pairing decides which combinations of the family's equations are kept. It matters for the distribution. The automatic pairing puts the mass equation among the equations of the distribution, and that equation is exactly the one a mass-preserving basis makes redundant. **Two stages.** ``'second_stage', true`` solves the reduced model of the first stage and compresses its coordinates by the decomposition above, family by family. Families kept exact are not touched; families with no basis, kept whole in the first stage, are compressed in the second. The bases compose, :math:`E = E_1 E_2`, and so do the projections, :math:`W' = W_2' W_1'`. By default :math:`W_2' = E_2' A_{0,1}^{-1}`, the projection above applied to the first-stage model at its own solution, which is exact on the span of :math:`E_2`. ``'second_stage_projection', 'galerkin'`` takes :math:`W_2 = E_2` instead. Under either projection the decomposition of the second stage includes the responses of every enforcer of the first-stage solution: :math:`W_2 = E_2` passes an enforcing shock to the reduced model only through the part of its response that the basis spans, the default projection is exact only on the span of :math:`E_2`, and the filters need that part to enforce the constraints. This is the pipeline of Bayer and Luetticke: a parameter-free first stage makes the model small enough to solve, and a decomposition of its solution compresses it further. On the model of the tutorial, the second stage takes the first stage's 103 coordinates to 49, the enforcers' responses included, with the first stage's likelihood. **Two tiers behind one interface.** The decomposition with ``'A0'`` adapts the basis to the solution: the reduced model approximates the full model, and it is exact whenever the full solution lies in the span of the basis. A supplied basis with ``'galerkin'`` needs no solution, so it makes a model feasible that could not be solved at full size. The reduced model is then an approximation of the full model when its basis spans the full solution, as with a complete basis. Otherwise it is a model in its own right, defined by the basis and the projection: compare its root count and its likelihood with those of richer bases. Both tiers produce the same object, and the filters, the constraints, regime switching and the higher orders apply to either. On the switching model of the tutorial, with constant and with time-varying transition probabilities alike: * a complete supplied basis reproduces the full likelihood to about :math:`10^{-13}` under every projection; * 8 DCT coefficients for consumption, and the marginals plus 5 copula perturbations for the distribution, keep the log-likelihood within :math:`4\times10^{-3}` of the full one with ``'A0'``, and within :math:`3\times10^{-2}` with ``'galerkin'``. Contracted derivatives ---------------------- At every solve, the reduced problem needs the derivatives of the model in the reduced coordinates, :math:`W^\top\, D_v\, E_v`, where :math:`D_v` is the Jacobian of the full model and :math:`E_v` maps the reduced differentiation variables to the full ones. A reduction built with ``build`` evaluates them **without assembling** :math:`D_v`: the rows of the symbolic derivatives (one value vector per template row under ``template_diff``, one entry per grid cell) are contracted straight into the reduced problem, in sparse chunks. The contraction is linear in the derivative values. The reduction keeps the values of the last evaluation, so a new one contracts only the entries whose value changed: .. math:: W^\top D_v^{new} E_v = W^\top D_v^{old} E_v + W^\top (D_v^{new} - D_v^{old}) E_v A draw that moves aggregate parameters only -- policy rules, shock processes, the switching probabilities -- changes the derivatives of a few equations and leaves the household block untouched: those rows are evaluated but not contracted. This is the device Bayer, Born and Luetticke build by hand (the household rows of their Jacobian are held fixed during estimation), here automatic and exact. .. code-block:: matlab R = rise.engine.dsge_tools.reduction.build(m); % 'contract', true m = set(m, 'solve_reduction', R); [~, LL] = filter(m); % first evaluation: every entry m = set(m, 'parameters', {'theta_pi', 2.2}); [~, LL] = filter(m); % only the entries theta_pi moves R.contract_cache.last % rows, entries, changed, full The result is that of the full evaluation followed by the projection, to rounding; the state is recontracted from scratch every ``R.contract_cache.refresh`` updates (200). On the 1,409-equation HANK above, the first evaluation gives the same log-likelihood as the full path to the last digit, with and without regime switching (constant and time-varying probabilities), and moving the monetary rule, the fiscal rule or the monetary shock contracts 1 to 6 of about 100,000 derivative entries. The contraction applies at first order with symbolic derivatives. Higher orders, the self-consistent linearization, and reductions built before the option existed use the full evaluation; ``build(m, ..., 'contract', false)`` opts out. Validating a reduction ---------------------- A reduction is an approximation. Before relying on it, compare the likelihood of the reduced and the full model at a few parameter vectors -- the reference, and points away from it in the directions an estimation will explore: .. code-block:: matlab for p = points [~, LLf] = filter(set(m, 'parameters', p)); [~, LLr] = filter(set(mr, 'parameters', p)); fprintf('%12.6f %+.2e\n', LLf, LLr - LLf); end At order 2 and above, compare **solutions**, not only likelihoods: simulate the full and the reduced model with the same shocks (``'simul_pruned', true``) and measure the gap against the higher-order part of the paths (the order-2 path minus the order-1 path). The nonlinear filters place their sigma points in the coordinates of the state they are given, so their likelihood changes with the coordinates: on the small model of the tutorial, by about 0.1 even with a complete basis (the full model written in other coordinates), while the simulated order-2 paths of the default reduction differ from the full ones by 6e-6 against a second-order part of 5. Two checks are exact and cost nothing to run on a small model: with a complete basis, ``build(m, 'tol', -Inf, 'forcing_tol', -Inf, 'enforcers', 'all')``, the reduced likelihood, impulse responses and smoothed variables equal the full ones to rounding error; and ``R.dimension`` reports the sizes (``full``, ``reduced``, ``kept``, ``coordinates``, ``states_full``, ``states_reduced``). On a HANK of 1,409 equations (the Bayer-Born-Luetticke household block on a 4x4x3 grid, 695 states), the default reduction has 485 variables and 162 states and reproduces the log-likelihood to within 0.03 at the posterior mean and at points 20% away in persistences, volatilities, policy rules and frictions, at 1.4-2.4 s per evaluation instead of 19-30 s. On a one-asset HANK with 100 borrowing limits that bind in sample, with and without regime switching, 228-242 variables keep the constrained log-likelihood within 0.02. Options of ``build`` -------------------- ``tol`` (``1e-10``) relative eigenvalue cutoff for the directions of the shocks; ``forcing_rho`` (``[0 .5 .8 .9 .95 .98 .99 .995 .999]``), ``forcing_periods`` (1500), ``forcing_tol`` (``1e-8``) the forcing of the kept block; ``enforcers`` (``'auto'``), ``enforcer_draws`` (5000), ``enforcer_scale`` (2) the enforcers of occasionally-binding constraints (see above); ``order_enrichment`` (``true``), ``second_order_directions`` (10) the higher-order enrichment; ``variables``, ``groups`` the variables to reduce and their families, when the defaults taken from the heterogeneity expansion do not fit; ``max_dim`` (``Inf``) the largest number of coordinates per group; ``basis``, ``projection`` (``'A0'``), ``equation_basis``, ``equation_families``, ``equation_matching``, ``second_stage`` (``false``), ``second_stage_projection`` (``'A0'``) supplied bases and projections (see above); ``contract`` (``true``) contracted first-order derivatives (see above); ``max_dense`` (``3e8``) the largest dense array of the reduced higher-order terms, in elements (see Higher orders). Limits ------ * Reduced variables must keep their steady state across regimes and parameter vectors; declared log variables, and variables with a balanced-growth component, cannot be reduced. * Observables are kept variables. * Conditional information inside the constrained filter (the conditional forecasting path) is refused. * The fast estimation kernel (``estim_fast_kernel``) solves the full system and declines when a reduction is set. * The filters report the means of the full variables; their covariances stay in the reduced coordinates. * A model object saved before the option existed cannot set it: rebuild it. See also :doc:`Heterogeneous agents` and :doc:`Very large models`.